clear;clc;clf;
p1=[];
p2=[];

k=500;
f=@(x)(k*x+1025)/(x^2+k);
for i=1:100
    x0=i;
    for j=1:100
        x0=f(x0);
        %fprintf('j= %i,x0= %f\n',j,x0);
        if x0^3==1025
           p1(i)=i;
           p2(i)=j;
            break;
        end
    end
end
plot(p2,p1,'r:.')

            